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DYNAMICS OF A ONE -DIMENSIONAL PLASMA SHEATH* 

By Frank Hohl and Leo D. Staton 
Langley Research Center 

SUMMARY 

Analytical and numerical methods are used to investigate a bounded one- 
dimensional plasma with a fixed neutralizing ion background. A variational method is 
used to show that a large class of stationary solutions of the nonlinear Vlasov equations 
represent minimum-energy states. The consequence of the minimum -energy property 
is that the stationary states, which represent minimum-energy configurations, can never 
be completely reached. However, numerical experiments which simulate the system 
reveal the interesting property that the plasma approaches its stationary state closely 
whenever the initial energy is not too different from the energy of this stationary state. 

INTRODUCTION 

A one -dimensional model is used to investigate the approach to an equilibrium 
state (that is, to a state described by a stationary solution of the Vlasov equation) of a 
bounded plasma with a fixed neutralizing ion background. The stationary states of the 
system are investigated analytically, and numerical experiments with a charge-sheet 
model are performed to study the time development of the system. The present problem 
is different from most previously investigated problems in that the system is strongly 
inhomogeneous. The fixed neutralizing ion background is constant over a given region 
and is zero outside this region. The electrons will then form a sheath at each side of the 
plasma slab. 

The usefulness of the one -dimensional sheet model has been well established in 
plasma physics. Some pioneering work in this field was done by Buneman (refs. 1 and 2). 
Buneman finds that collective interactions or collisions of charge in bulk can restore a 

*The information contained in the text is part of a thesis by Frank Hohl entitled 
"Collective Effects in Stellar Dynamics and Plasma Physics" offered in partial fulfill- 
ment of the requirements for the degree of Doctor of Philosophy in Physics, College of 
William and Mary, Williamsburg, Virginia, June 1967. The appendix was written by 
Leo D. Staton. 



grossly non-Maxwellian velocity distribution within a few plasma periods to a near- 
Maxwellian distribution by means of the buildup of electrodynamic instabilities. In 
similar work, the one -dimensional sheet model has been applied to specific problems 
by Dunn and Ho (ref. 3), Birdsall and Bridges (refs. 4 and 5), Smith and Dawson (ref. 6), 
Burger, Dunn, and Halsted (ref. 7), Derfler (ref. 8), and Hasegawa and Birdsall (ref. 9). 
The same model was used by Burger (refs. 10 and 11) to analyze the stability of direct- 
current solutions for the plane diode. A discrete one -dimensional computer model for a 
collisionless plasma in a magnetic field was used by Auer, Hurwitz, and Kilb (refs. 12 
and 13) to examine the structure of magnetic compression waves. This work was later 
extended by Rossow (ref. 14). 

Dawson (refs. 15 and 16), and Eldridge and Feix (refs. 17 and 18) initiated a study 
of the one -dimensional sheet model to check the theoretical predictions of plasma 
behavior. Dawson (ref. 15) investigated the thermal equilibrium properties of plasmas, 
such as: the drag on a test particle, Debye shielding, diffusion in velocity space, Landau 
damping of Fourier modes, and other parameters. In a subsequent paper (ref. 16), 
Dawson used numerical methods to investigate the thermal relaxation of a one-species 
one -dimensional plasma and found that there is no relaxation to a Maxwellian to first 
order (that is, a number of plasma periods equal to the number of particles per Debye 
length), as was shown analytically by Eldridge and Feix (ref. 19). Eldridge and Feix 
(refs. 17 and 18) performed numerical experiments to study the properties of the one- 
dimensional plasma near equilibrium and related some of the one -dimensional plasma 
properties to the three-dimensional case. The agreement with the theoretical prediction 
of the results obtained from the numerical experiments gives considerable confidence in 
the simulation of Vlasov systems by only a few thousand particles. 

SYMBOLS 


a acceleration 

a nm matrix elements defined by equation (79) 

N 

A constant, - 

4x 0 v 0 

B constant defined by equation (64) 

e magnitude of electronic charge or charge per unit area 

E electric field 


2 


f 


electron distribution functions 


F energy distribution function 

g energy density defined by equations (72) and (76) 

g* = g + 2AB k v( k ) 

g-j> kinetic -energy density, B k v( k )^ 

h function defined by equation (85) 

H Heaviside unit step function 

m electron mass or mass per unit area 

n density 

N total number of electron charge sheets 

P potential energy 

s = 1 - z 
t time 

T kinetic energy 

1 9 

U total energy of a charge sheet, -mv - ecp 

Ct 

v velocity 

v-j, thermal velocity 

V,V contours defining regions of constant f 

w variable defined by equation (A4b) 

W total system energy, T + P e + Pj 
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w* 

x 

x,x 1 ,x 2 

z 

a 

P 

T 

8 

e 

e f 

C 

V 

e 

K 

A 

X D 

P 

a 

T 


defined by equation (A5) 


position coordinate 


positions defined by V(x) = 0 


dimensionless variable, 


Nef 

8ef€ 


constant defined by equation (25) 


2A / 7 r 

constant of integration, ^on 

defined by equation (A4c) 


Dirac delta function or dimensionless energy, y 


maximum energy of an electron defined by equation (38) 


permittivity of free space 


dimensionless position coordinate defined by equations (22) and (49) 


dummy integration variable 


variable defined by equation (70) 


thermal energy 


Lagrangian multiplier 


Debye length, v-p 


(me. 


n e e 


charge density 

charge per unit area 

period defined by equation (45) 
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<p 


electric potential 


<3? dimensionless potential defined by equations (21) and (50) 


n e^ 

top plasma frequency, — — 


Subscripts: 


e 

electron 

eq 

equilibrium 

i 

ion 

j>k,n 

summation indices 

min 

minimum 

o 

initial 

s 

system boundary 


Superscripts: 

k summation index defining waterbag contours or ordering index 

m,n summation indices 

A prime on a symbol indicates partial differentiation with respect to x. 

STATIONARY EQUATIONS FOR A BOUNDED PLASMA 
The characteristic time for collective effects in a plasma is given by 



where e is the charge and m is the mass of an electron. The characteristic length 
is the Debye length given by 


I 



(lb) 



The Vlasov equation for the one-dimensional system with a fixed ion background 
takes the form 


9f + v 9f_ _ eE 9f _ 0 
at 9x m a v 


( 2 ) 


where E = -^ and f is the electron distribution function. The electric potential cp 
9x 

is given by the Poisson equation 



(3) 


where nj is the fixed ion density. The integrations are performed over all values of 
the variables throughout the remainder of the paper, unless otherwise specified. In the 

r)f 

steady state, — = 0 and equation (2) takes the form 

V !1.£E 9f_ =v _9_ + _e d£8L =0 ( 4 ) 

9x m a v 9x m dx 9v 


Virial Theorem 

In order to obtain a relation between the kinetic and potential energy of the system 
in equilibrium the identity 


-2L. y ^ x2 f(x,v)dx dv = m ^ ^ v^ f(x,v)dx dv + m ^ ^ xa f(x,v)dx dv 


(5) 


is investigated where 


D 9 9 9 

— = — (. v f a — 

Dt 9t 9x 9v 


( 6 ) 


and — = 0 according to equation (2). Because the system is in equilibrium, the left - 
Dt 

hand side of equation (5) is zero. The first term on the right-hand side of equation (5) 
is 2T where T is the expectation value of the kinetic energy of the electron system. 
The last term of equation (5) is now investigated. By using Newton's law, a can be 
replaced by ^ Also, the electron density is given by 
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n e = J f dv 


(7) 



By using the Poisson equation 


dx 2 




the electron density can be written as 


n e 


6 f d 2 cp 
e dx 2 


+ m 


( 8 ) 

(9) 


The last integral of equation (5) can then be written in the form 



An integration by parts gives the result 




( 12 ) 


where P e is the potential energy of the electron system and the system is assumed to be 


symmetric in x. The term 


in the potential-energy definition is required 


-x c 


to normalize the potential energy if the system contains a net charge or if an externally 
applied field is present. The last integral in equation (10) must now be considered. An 
integration by parts gives the result 


r x i 

\ en i: 




x dx = en^xd 
dx 1 


p x i 

- \ en^cp dx = -Pj 


(13) 


-X< 


where Xj defines the boundary of the fixed neutralizing ion background; that is, 
nj = n G for -xj < x < Xi and n^ = 0 otherwise. The energy P^ is related to the 
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potential energy of the ion background. In general, the potential energy of a one- 
dimensional system of charges 0 } is given by 


p = i £ m - f y PC’ 


( 14 ) 


Because <p = 0 can always be chosen at x = xj, equation (13) can be put in the same 
form as equation (14). In the remainder of the paper cp is always taken to be zero at 


x = 0. Only systems with zero net charge are treated so that x 
tion (6) takes the form 




-x s 


= 0, and equa- 


2T = P e + Pi (15a) 

In the limit as X| — 0, the system consists of two electron sheaths which are held 
together by a thin but dense positive charge. The ions then simply take the place of a 
positively charged grid through which the electrons can pass freely and Pi = 0. Equa- 
tion (15a) then takes the usual form 


2T = P e (15b) 

The total energy of the system is always given by T + P e and is a constant of the 
motion. 

If the ions are not assumed to be fixed and the ion density is given liy an equation 
similar to equation (7), then 


2(T e + T t ) - P = 0 (16) 

where P = P e as given by equation (12). 

Equilibrium Solutions 

By using the method of characteristics to solve equation (4), any distribution func- 
tion F, which is a function only of U, is found to be a solution of equation (4) where 

U = i mv 2 - e<p (17) 

Ci 

is the energy of a charge sheet. If an equilibrium state can be reached by starting from 
a given initial distribution f 0 , then 
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f 0 (x,v,t - °°) - F(U) 


(18) 


The form of F(U) depends on the initial distribution and must, in general, be obtained 
by following the time development of the nonlinear Vlasov equation. However, if an 
initial distribution which is constant over a certain region of phase space and is zero 
outside this region is taken, then F(U) can be obtained without following the evolution 
of the distribution function. This F(U) is unique. The distribution is shown in figure 1. 
According to the Vlasov equation, f remains constant along the different trajectories 
so that the region in phase space can only change its shape with time while its area 
remains constant. For this reason, this distribution function has been called the water - 
bag model by DePackh (ref. 20). This model is of interest primarily because there 
corresponds to it only one possible equilibrium state which is easily calculated and can 
be used for comparison with the numerical results. 

For the waterbag model, F(U) = A = Constant for 0 SU se, and F(U) = 0 
for U > e . 

The equilibrium solution for the waterbag model is now obtained. The initial 
shape of the waterbag is taken to be a rectangle defined by the area in phase space 
between ±x 0 and ±v 0 . The function F(U) is used in the Poisson equation (3) and 
the integration over dv is changed to an integration over dU. Because the density at 
a given x is of interest, dv is given by 


±dv = 


dU 

\! 2m(U + ecp) 


(19) 


The Poisson equation can be written as 



where 

n i= n o (|x| < |Xi|) 

and 


n^ = 0 


(l x l > l x il) 
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If 4> and C are introduced such that 



( 21 ) 


and 



then equation (20) takes the form 

- 1 + a.\) 1 - 4 > = 0 

d? 2 

and 

^-^■+ ot\J 1-4? = 0 

d? 2 

where 


Q! = 


2A fgT 
n Q Vm 


A first integration of equations (23) and (24) yields 


( 22 ) 

(? < ?i) (23) 

(C > ?i) (24) 

(25) 


and 



(?<q) ( 26 ) 

(e > ?i) (27) 


respectively. The constants of integration in equations (26) and (27) have been evaluated 
by using the boundary conditions = 0 and $ = 0 at C = 0, and = 0 and 4? = 1 

at C = Co where C defines the boundary of the system. The first two boundary con- 
ditions are obtained by applying symmetry in C an d the last two are obtained by using 

J T 

overall charge neutrality. Because — and <3? are continuous at the boundary of the 

dC 

neutralizing ion background Ci 

$(C i ) = $ i = ^ (28) 
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II 


and a second integration gives the final result 



where ^ is given by equation (29) with the upper limit of the integral equal to <3?^. 

The numerical solution of equations (29) and (30) is shown in figure 2 which gives $ as 
a function of In £ for several values of a. Equations (21) and (22) can be used to 
obtain <p as a function of x. The relation 


n 0 = 


N 

2xi 


(31) 


where N is the total number of electrons or ions in the system and equation (25) is 
used to eliminate n 0 and e from equation (22). Then 


x i = 




(32) 


and the numerical values for n Q and e are easily obtained. The constant of propor- 
tionality in equation (22) is found from the ratio 


ee. 


n o e 


x 


_2 

Si 


(33) 


where equation (32) is used. Figure 3 shows the variation of the electron density for 
three values of a for a system with N = 1000, A = 0.25, and = m = e = 1. The 
variation of the electric field is shown in figure 4 for corresponding values. The 
equilibrium contour in phase space V ± (x) of the waterbag is given by the equation 

V + (x) = -V_(x) = J(e + ecp)^ (34) 


and figure 5 shows the equilibrium contours for three values of a. The parameters 
used in the calculation are the same as those used in figure 4. Equations (29) and (30) 
cannot be used if the positive ion background reduces to a thin but dense charge and 
Xj — 0. Such a system corresponds to a positively charged permeable grid through which 
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the electrons can pass freely; that is, there are two electron sheaths that are held 
together by a dense, thin positive charge. The Poisson equation (3) can be written as 

(x/0) (35) 

dx z 6jv2m 

A first integration gives the result 


d .(£_ _ 
dx 


16A : (e + ecp + Constant 


f 3e f /2m 


(36) 


The boundary conditions are that ^ = 0 and cp = - — for x = x~ and that 

dx e ° 

lim ^ and lim cp = 0. By using these boundary conditions, the constant 

x-0+ dx 2 e f x— 0 + 

of integration is found to be zero and 


d£ = _Ne_(i e£\ 
dx 2e# \ e / 


3/4 


(37) 


where 


; 3/4 


/ 16A Ne 

3e f /2m 2e 


f 


so that the maximum energy of an electron is 


An integration of equation (37) gives the final result 


-\4/3 


(38) 


dx = 


8ej£ 

Ne 2 




1/4 


(39) 


Figure 6 shows the variation of E = - ^ given by equation (37) as a function of x. 

The variation of the period of oscillation of a charge sheet as a function of energy 
seems to give a measure of the effectiveness of phase (or orbit) mixing for the system. 
For example, in the problem of a self -gravitating system (ref. 21), the variation of the 
period with particle energy is very small and the oscillations in the kinetic energy of the 
system persist for a long time. As is shown subsequently, the variation of the electron 
period with energy is large for the present problem and the oscillations of the kinetic 
energy with time quickly damp. Only the period for the case where Xj = 0 will be 
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calculated. The equation of motion for a charge sheet is 


d 2 x _ e d<£ _ _ e d 
2 dx dx 


( 40 ) 


By using the normalization given by equation (21), equation (39) (for positive x) can be 
written as 


# = 1 - (1 - z) 


(41) 


where z = x 


Ne 

,8e f e 


. Therefore, 


4(1 - z) 3 
dz 


(42) 


and equation (40) takes the form 

<^ + _NV (1 ..)3 = 0 

dt 2 16m€ f 2 e 

A first integration of equation (43) gives the result 


/dz\ 2 _ N 2 e 4 [, 

i ji I o 1' 

1 - z) 4 - (1 - 6)1 

/ 32me f 2 e L 

— 1 


(43) 


(44) 


where 6 = and U is the total energy of the sheet under consideration. A second 


integration of equation (43) gives the period of oscillation as 


16e 


r = 


Ne 

16e 

Ne : 


n l-(l-5) [- . 

2 \/2me j |_(1 - z) 4 - (1 - 6)J dz 

■ \[2me C 1 Is 4 - (1 - 6)1 ^ ds 

(1-6) V 4 


(45) 


where s = 1 - z and the limits of integration are found from equation (44) by setting 

— = 0. The variation of the electron period as a function of energy is shown in figure 7. 
dt 

The time t q has been arbitrarily chosen. Because of the strongly nonlinear restoring 
force eE(x), the variation of the period from zero to infinity is not surprising. For the 
case where xj is not zero, the period approaches infinity as U approaches e . 
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The stationary state for a system which is described by a distribution function 

F(U) = A exp(-xU) (46) 

is easily obtained by the same method as that used for the waterbag model. The solution 
is given by 


and 


f r -I 

, -1/2 



5 = J o |2 [t? + /3 exp(-? 7 )J - 2/3 

f d?7 

(? < k) 

(47) 

? = yi" exp (f-) - exp (ij 

+ ?i 

6 > <0 

(48) 


where 


J8 = 


2A rir- 
n o V 2/cm 


■. Also, 


and 


2 = e 2 n 0 x 
e f 


(49) 


§ = -ecpK 


(50) 


were used. All but one of the constants of integration can be directly evaluated by using 

the boundary conditions ^ = 0 and <f> = 0 at £ = 0, lim — = 0, and continuity of 

d ? £-oo d? 


— at £ = The equation for £ > q 
constant is nonzero. In order to satisfy 


is then integrated by assuming that the 

j X 

lim — — = 0, the integration constant must be 

£-00 d£ 


zero, and equation (48) is obtained. For the case xj = 0, the solution is 



(51) 


(52) 


The constants of integration have been evaluated by using the boundary condition 

lim = 0, lim + ^ and <3? = 0 at x = 0. The condition lim — = ^ e2/c 

x-oo dx x— 0 dx 2e f x-0 dx 2e f 


is 
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obtained by applying Gauss’ theorem to the positive charge sheets at the origin. The 
electron density corresponding to equation (51) is 


n e (x) 


N 2 e 2 /< 


8e 


f 



( 53 ) 


Figure 8 shows $(x) and E(x) for the case where Xj = 0, e f = 1, k = 1, e = 1, and 
N = 4. The electric field is very similar to that for the waterbag case. Thus, irrespec- 
tive of the exact form of the distribution function, similar evolution of the system might 
be expected. 


MINIMUM-ENERGY PRINCIPLE 


The stationary solution for the initial waterbag distribution was obtained by con- 
serving only area in phase space. The stationary (or equilibrium) state is described by 
a function which depends only on the energy U. For the waterbag model, a unique F(U) 
corresponding to a given initial energy and particle number can be obtained. Before 
such an F(U) can be reached, energy considerations must allow the relaxation from a 
function f(x,v,t = 0) of two independent variables to a function F(U) of the energy 
only. If an equilibrium state is to be reached, the initial and final energy of the system 
must be equal. For the initial rectangular waterbag, the kinetic energy T is given by 


T(t = 0) = J J | mv 2 f(x,v,t = 0)dx dv 

= Nm_ f X ° f V ° v 2 dv ^ 

8x 0 v 0 *J-x 0 J-v 0 

Nm 2 
6 ° 


(54) 


As previously stated, ±v 0 and dx 0 define the rectangular area of the initial distribu- 
tion. For positive x, the electric field is given by 


and 


E(x) = 


Ne / 1 

2e f \ x i 



(x < Xi) (55) 



(*i < x < x Q ) (56) 
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and, therefore, the initial potential energy is 


p e (t 



E 2 (x)dx 


e 2 N 2 

12ejX 0 




( 57 ) 


Again, ±x Q and ±v Q define the rectangular area of the initial distribution. The total 
energy of the initial distribution is given by 


W(t = 0) = T(t = 0) + P e (t = 0) 


N 

6 


mv 2 + -Nfii 

° 2 Vf 


( x o - x i)‘ 


(58) 


The minimum value that equation (58) can attain for a given value of A = 

by 


N 


4x o v o 


W 


mm 


(t = 0) = 


Nm j 


( N 2 e 2 
1 16 Aejm 


2/3 


is given 


(59) 


where x^ = 0. 

By using equations (37) and (38), the potential energy of the equilibrium state for 
xj = 0 is given by 



16e 2 / Ae f 
" 7e y 3 /2m7 


The kinetic energy is given by 

px s pV + (x) x 

T en = 4 \ \ imv 2 A dv dx 

eq Jo Jo 2 

_ 8e 2 / Ae f 

V 3 \j 2me 


(60) 


(61) 
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and the initial rectangular distribution with the smallest possible energy has still more 
energy than the equilibrium distribution. The equilibrium waterbag distribution for a 
plasma sheath may have, therefore, a minimum-energy, property. 

The minimum -energy property for the stationary solutions of a plasma sheath 
represented by a single waterbag is given in the appendix. This property can be shown to 
hold for a large class of distribution functions that can be described by multiple-contour 
waterbag distributions. The waterbag model illustrated in figure 9 is used in the analy- 
sis. The contours V + ( k )(x,t) and V_( k )(x,t) describe surfaces of constant f = fj^. 
According to the Vlasov equation, f behaves like an incompressible phase fluid; there- 
fore, the area bounded by any pair of contours is constant. With a large number of con- 
tours, the waterbag model can be used to approximate arbitrary distribution functions. 

The function 

f = ^ B k |h [v - V_( k 5 - H [v - V + ( k > | (64) 

k 

describes the waterbag distribution where the summation is over all contours and H(z) 
is the Heaviside unit step function. The distribution function f, as given by equation (64), 
satisfies the Vlasov equation 
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where 6(z) is the Dirac delta function. The electric field E is obtained from the 
equation 


dE 

dx 




( 66 ) 


where there is assumed to be a total of N electrons each of charge e in the system. 
If equation (65) is integrated separately over positive and negative velocity 


and 


I 


B,. 


8 V 


(k) 


at 


■+ V 


(k)!X± 


(k) 


8x 


eE | 
m 


= 0 


l 


9V_ (k) ( k ) 8V_ (k) [ eE 


at 


• + v 5 


8x 


m 


= 0 


(67a) 


(67b) 


The stationary contours are described by 



k 


(68) 


which must hold term-by-term. In order to simplify the equations, symmetric contours 

V + ( k ) = _v_( k ^ = V k ^ are assumed in the remainder of the analysis. The results 
obtained are not affected by such an assumption. Equation (68) can then be simplified as 


y (k) dV^. + eE = Q 
dx m 


(69) 


A new variable is now defined as 


6>( k \x) = 0 


(*<«.«) 


n(k) 


(x) = r V (k) (C)d? (-x s (k) Sxl x s ^) 

J (k) 


-x, 


> ( 70 ) 


» ,k >(x) = 8 <k) [x s < k j] 


( x >x s W) 


J 
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where ±x s ^) are the end points of the kth contour. In terms of 0^) the electric 
field is given by the equation 



H^x + XjJ 


H 


(x - x^jj + X| H^x + x^ - H^x 



( 71 ) 


The energy density of the system can be represented by a function 
total energy of the system is then given by 



The 


W = 



(72) 


For the integral for W to be an extremum subject to the condition that 


N = 


I 


2B k 0^ 



(73) 


is a constant requires that the contours 



satisfy the Euler -Lagrange equation 

d._^L _ =0 

ay( k ) 


(74) 


where g* = g + X2B k V^ and X is an undetermined Lagrangian multiplier. Because 

the end points ±x s ^ are fixed, the constraint on the number of particles is redundant 

and X = 0. If the end points ±x s ^ are not held fixed in the variational process, addi- 
tional equations which occur do not change the main results. 

For a multiple waterbag, the kinetic energy per unit length is given by 


X g T & (k !] = I f B k y(k)3 ( 75 > 

k k 
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The expression for g then takes the form 

,3 e. 


l g[x,V«,^ - l & V « 3 * {x[H(x . x t) - H(x - Xj) ] 

k k 

+ Xi [H(x + Xj) - H(x - Xi )J j - Y B k e (] 

J k 

= J ^-V (k)3 - 2 e £ B k e^h^- fx [H(x + Xi ) - H(x - x^ 
k k ' 

+ Xl [H(x + Xi) - H(x - Xi)] j - Y B k e (k) \ 

+ t(?) 2 { x [ h ( x + Xi ) - H ( x ‘ x i3 

2 

+ Xi |h(x + x^ - H(x - Xj) \ ( 76 a 

tk) 

Because the last term in equation ( 76 a) does not depend on the function V v ’ or 
8^\ it can be omitted and g is thus given by 

£ g[^,V (k) , 0 (k j] - £ 5|^ V ( k)3 - 2 e |x|H(x + Xi ) - H(x - Xi )J 


+ x i [ H ( x + x i) - H (x - Xi)] j - ± £ B k 0 (k) j 


(76b) 


The Euler -Lagrange equations then become 


2B l 


mV 


« ^ + eE 
dx 


= 0 


(77) 


or 


y (k) dV (k) + eE = Q 
dx m 


( 78 ) 


which is the equation for the equilibrium contours. 
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If equation (78) is to represent a minimum -energy configuration, then Legendre's 
criterion of the second variation of g must be satisfied. The Legendre condition 
(ref. 22) for the case of several unknown functions 0^) j s that the quadratic form the 
coefficient matrix of which has the elements 

a = (7« 

11111 9y( n )av( m ) 

must not be negative. For the present problem, only the diagonal elements of equa- 
tion (79) are nonzero and are given by 


a^ = 2mB k V^ (80) 

Because 2mV^ is never negative, the Legendre condition requires that 

B k > 0 (81) 

The form of f given by equation (64) shows that the condition >0 is equivalent to 
stating that the distribution function must decrease monotonically in going outward from 
the center of the system where f = f^ = f^ must be the largest. If equation (81) is sat- 
isfied, the system is not a maximum -energy configuration but is a minimum -energy 
configuration which is therefore stable. However, if > 0 is not satisfied for all k, 
the system is not a minimum-energy configuration and nothing can be said about its 
stability. 

If the distribution of ions is not fixed but is given by an f similar to equation (64), 
the system is found to be a minimum -energy configuration when both the electron and ion 
distribution functions are monotonically decreaseing in going outward from the center of 
the system. 

For the case of a single -contour waterbag, the Legendre criterion is always satis- 
fied and the stationary solution represents a minimum -energy configuration which is 
inaccessible for arbitrary initial distribution. However, the results of numerical exper- 
iments show that the system approaches the equilibrium state closely whenever the ini- 
tial energy is not too far from the energy of the corresponding stationary state. The 
single -contour waterbag is treated in detail in the appendix. 

THE ONE -DIMENSIONAL PLASMA MODEL 

The plasma model consists of a system of N ions and N electrons that are 
represented by 2N charge sheets. The charge sheets move only in the x-direction and 
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are of infinite extent in two transverse directions. When two sheets meet, they are 
allowed to pass freely through each other. The equations of motion of all the 2N sheets 
are solved simultaneously by computing the position of each sheet and by integrating the 
equations of motion of each sheet over a small time interval At. The equation of motion 
of a sheet with charge Oj and mass mj per unit area is given by 


d^^ 
dt 2 m j 3 


(82) 


where Xj is the position of the sheet. The electric field is given by 

9E(x,t) . 1 , t) 

9x e f Mx ’ g 


2N 

= ^6[x-xi(tj] 


(83) 


i=l 


The magnetic field induced by the charge sheets is neglected. The electric field is then 
given by 


where 


E(x -‘ ) = h Z°‘ h [ x ■ Xi(t il 


(84) 


h(z) = 1 
h(z) = 0 
h(z) = -1 


(z > 0) 
(z = 0) 
(z < 0) 


(85) 


In obtaining equation (84), the electric field is assumed to be zero outside the system. 
For the numerical calculations, the electric field acting on a particular sheet is obtained 
by the net charge to the left of the sheet Xj; thus 



_1_ 

e f 





( 86 ) 
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where the superscript k gives the order of the charge sheet xy ' such that 


x (k) < x (k+l) 


(87) 


The IBM 7094 electronic data processing system was used to calculate the self- 
consistent motion of systems containing several thousand charge sheets. The position 
and velocity of each sheet are computed at successive times tj, tg, tg, . . . t m . For 
each time step At = t n+ i - t n , the new positions and velocities of the sheets are com- 
puted from the equations 




( 88 ) 


and 


_ <fri(*n) , d2x i( t n) At + 1^W (M ) 2 

dt dt dt 3 


2 dt 3 


(89) 


where 


and 


°i 

dt 2 e f m i 


I °' (k) - k a i 

k=l 


(90) 


dt c 


1 

d 2x i(t n ) 

d2x i(Vl) 

At 

dt 2 

1 

CM 

-4-> 

1 


(91) 


The kinetic energy of the system is given by 


y m. /dx- N 2 

L 2 \ dt 

3=1 


(92) 


and the potential energy is 


2N-1 


p ’i I (Vi- x i) E j : 

3=1 


( 93 ) 
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NUMERICAL RESULTS 

All numerical computations for the sharply bounded plasma were performed for 
the case where the ratio of the ion mass to electron mass is so large that the ions can be 
assumed fixed and only the motion of the electrons needs to be calculated. The problem 
is that of determining the evolution of the electron distribution function for a fixed ion 
background between -Xj and x^. The case where xj is near zero is investigated 
first. The case corresponds to that of a permeable grid which is positively charged and 
through which the electrons can pass freely; that is, there are two electron sheaths held 
together by a thin but dense positive charge. For all calculations, e=ej = m= l is 
chosen. Figure 10 shows the variation of kinetic energy for a 1000-electron system. 

The ratio of the initial energy to the equilibrium energy, as given by equation (71), is 1.3. 
The dashed line shown in figure 10 corresponds to one -third of the total energy of the 
system and is the kinetic energy given by the virial theorem. The kinetic energy quickly 
approaches this value. The plasma frequency o>p 0 corresponds to the initial electron 
density. One important indicator of the approach to an equilibrium state is the time 
development of the energy distribution function. Figure 11 shows that, after about six 
plasma periods ^ 27 ra>p 0 ~-^ , the energy distribution function has approached the equilib- 
rium distribution indicated by the dashed line. After six plasma periods, the system 
changes very little as its evolution is followed for a longer time. The high-energy tail 
of the distribution does not disappear and is required to conserve the energy of the sys- 
tem. The corresponding time development of the density is shown in figure 12. The 
time development of the density indicates that the density of the system quickly tends to 
approach its equilibrium value as indicated by the continuous curve. The dimension of 
the sheath at each side of the system is of the order of a Debye length so that the number 
of particles per Debye length is of the order of the number of particles in the system. 
Dawson (ref. 16) has shown that thermalization in a uniform one -dimensional plasma 
occurs only on a time scale Graininess effects in the present system can 

be expected to become important only after 10® plasma periods. The results shown in 
figure 12 are for times much less than (nAp^Wp and are therefore not affected by 
graininess. The electric field at t = 30u>p O "l is compared with the equilibrium field 
in figure 13. Because of the high-energy electrons, the electric field of the system does 
not go to zero as quickly as the equilibrium value. The time development of the system 
in phase space is shown in figure 14. The system quickly develops arms which wind 
around the main body of the system as it rotates in phase space. Because of the very 
long period of rotation for particles at the outer boundary of the system, the arms 
stretch rather quickly. The calculations in figure 14 are repeated for a 2000-electron 
system. The results are almost identical and only the initial time development of the 
system in phase space is shown in figure 15. 
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Calculation for values of x^ ^ 0 are also performed. The dimensions of the 
sheaths which form at each side of the positive ion background are of the order of a 
Debye length. For appreciable values of x*, most of the electrons are in the main body 
of the plasma, and the number of electrons in the sheath or per Debye length becomes 
quite small. For moderate values of xj, graininess effects become very important and 
the system approaches a Maxwellian distribution. Figure 16 shows the variation of the 
kinetic energy for a 1000-electron system with xj^ = 477. Figure 17 shows the time 
development of the energy distribution function for a collision-dominated system and 
shows that the system quickly approaches the Maxwellian distribution as indicated by the 
dashed-line curve. The corresponding results for the time development of the system 
in phase space are shown in figure 18. Figure 18 indicates that at t = 45.00wp O ~l the 
Vlasov character of the system has been lost and the system cannot be considered colli- 
sionless; therefore, the evolution of the system cannot be approximated by means of the 
Vlasov equation. 


CONCLUDING REMARKS 

The stationary state for the single -waterbag distribution is shown to be a minimum - 
energy configuration for the one-dimensional bounded plasma with a fixed neutralizing 
ion background. Thus, the stationary state is not accessible for any system which is 
initially in a nonequilibrium state because the initial distribution has an excess of energy 
which cannot be accommodated by the stationary state. Numerical experiments with a 
one-dimensional charge-sheet model reveal the interesting property that the system 
does its best within the limitations of energy conservation to approach the steady state. 
Similar results were previously found for a one -dimensional self-gravitating system. 

The minimum -energy property was extended by approximating arbitrary distribu- 
tion functions by a multiple -waterbag distribution. Stationary distribution functions, 
which decrease monotonically in going outward from the center of the system, are found 
to be minimum -energy configurations and are therefore stable. The consequence of the 
minimum -energy property is that for arbitrary initial distributions, described by a mul- 
tiple waterbag, an equilibrium state (in the Vlasov sense) is in general, inaccessible. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., January 24, 1968, 

129-02-01-01-23. 
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APPENDIX 


MINIMUM -ENERGY PRINCIPLE FOR A SINGLE WATERBAG 


In order to discuss in detail the minimum -energy property of a single-waterbag 
equilibrium distribution, the extremum of integrals of the form 



must be considered subject to the condition that 


(AI) 


N = 



V dx 


(A2) 


is a constant. The energy density g is a known function of X, V, and 0, where 
r X 

0 = j V(x)dx, and Xj and Xg are not fixed at the outset but V(x) must vanish at 

the end points. The positive constant A is equal to the height of the waterbag distribu- 
tion function. 

The purpose of this appendix is to establish the conditions that the smooth function 
V(x), which describes the equilibrium waterbag contour, must satisfy in order that the 
quantity W, which is subsequently identified with the total energy of the system, shall be 
an extremum or, in particular, a minimum. For this purpose, the one parameter family 
of integrals 


pX 2(6 ) 

W(e) = \ g(x,w,T)dx (A3) 

J X 1 (e) 


and 


N(e) = 2A 


P x 2 (e) 

^X x (e) 


w dx 


(A4a) 


are constructed, where 


W = V(x) + 6 7j(x) 


(A4b) 
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and 


T = 0(x) + e ?(x) (A4c) 

The function V(x) (as yet unknown) is the actual extremal curve and passes through 
the (as yet unknown) actual end points x^ and X2; t?(x) is an arbitrary differentiable 
function; e is the parameter; and £(x) is the indefinite integral of The trial 

function w(x) approaches V(x) as e approaches zero. The intersections of a 
given w(x) with the x-axis define the points and X2. As w(x) approaches 

V(x), the points Xj and X 2 approach the actual end points x^ and x 2 , respectively. 

By the usual arguments of the calculus of variations (ref. 23), the value of W(e) 
in equation (A3) for a given rj(x) is an extremum when e is zero. In order to take the 
constraint equation (A4a) into proper account, the integral 

pXo(c) 

W (e) = \ g*(x,w,r)dx (A5) 

^X 1 (e) 

is constructed, where g* = g + A(2Aw) and A is an as yet undetermined Lagrangian 
multiplier. The condition for an extremum is then 

= 0 (A6) 

S£ e=o 


With the use of equation (A5) 


9W* 

9e 


^g*[x 2 ,w(x 2 ), r(x 2 ] 



8^8w + 8 s l8T\ dx 
9w de dr de J 


(A7) 


Also, — = 77 and — = £ so that by using these relations in equation (A7) and by 
9e 

letting e = 0 

9W* 

de 

+ f ^ rj + Adx = 0 (A8) 

-’x. \ 9V 80 I 


6=0 


9X 2 

de 


g*[x 2 ,V( X2 ),«(x 2 )J - s. g*[x 1 ,V(x 1 ),0(x i y 
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8X 2 

The end-point conditions must be examined in detail to evaluate — - — 

9e 

9X! 

9e | e= o 


e=0 


and 


Because w(XjJ = w(X 2 ) = 0 by definition, then 

9w ( X l) _ Q 
9e 

for all e . Thus, the use of equations (A4b) and (A4c) gives 

9w(Xi ) q-k, 9Xi . . 9Xi . . , . 

v '( X l)-5T + ”( x l) + e " ( x l) = 0 

For e = 0 


(A9) 


(A10) 


V 


’( X l) 


9Xi 


9e 


6=0 


H X l) = 


= 0 


or 


By similar reasoning 


axi 77( Xl ) 

e=0 V( Xl ) 


3X 2 = _ Vfe) 

96 e=0 V '( x 2) 


(All) 


(A12) 


Substitution of equations (All) and (A12) into equation (A8) and use of integration by parts 
in the first integral of equation (A8) gives 


g*[ X 2, v ( x 2),e(x 2 )] g*|*l> v ( x l)> e ( x l)] 



v '( x 2) 

+ r z 

9 g * _ d fag\ 

4, 

90 dx \ 9V / 


^( x 2) 


£ dx = 0 


V 


'( X l) 


r?( x l) + ? 


9V 


(A13) 


Equation (A13) must hold for a rather wide class of functions rj(x) and £(x). In 
particular, equation (A13) must hold for those functions C(x) which are nonzero but are 
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compatible with = ?( x 2 ) = ^( x l) = ^( x 2 ) = °* For such a case 


Ml . i_(Ml\ = o 

90 dx\9V / 


(A14) 


Again, 77 can be nonzero at one end point but zero at the other end point when £ van- 
ishes at both end points, and so forth. By this method 


g*[x 2 ,V( x 2 ),e( x 2)] _ Q 

V t (x 2 ) 


(A15) 


g*[xi, v (xi),0(xi)] 

V>(xi) 


(A16) 


Ml 

9V 


= 0 


X=Xi 


and 



9V 


x=x 2 


= 0 


(A17) 


(A18) 


Equations (A 14) to (A18) and the original constraint equation (A2) with x^ and x 2 
replacing Xj and X 2 , respectively, are the equations used to determine V(x), xj., 
x 2 , and A. 

The function g(x,V, 0) for specializing to the case of the single waterbag is now 
introduced. The total energy per unit length of the waterbag g(x,V, 0) is given by 


f ±mv 2 f dv + e -^-= f V(X) A(imv 2 W + ^ E 2 

J 2 2 J_ v(x) \2 ) 2 

= ~ ^ jn Q Jh^x + Xi^ - H(x - x^jj - 2AV(x)J dx 

= Tr{ X [ H ( X + Xi ) " H ( x " x i)] +x i[H( x + x i) " H ( x - x i)]} “f ^ A0 


where 


(A19) 
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Then 


g*(x,V, 9) = g + 2AXV = V 3 + y(~T^ {* H ( x + x i) - H ( v * * * * x - x i)J 
+ x^|h(x + x^ - H^x - Xj^ J ^ A0^ + 2AAV 


(A20) 


Use of equation (A20) in equations (A17) and (A18) gives 


X = 0 


(A21) 


where V^x-jJ = V( x 2 ) = 0* For finite V'^x^^ and V'^x 2 ^, equations (A15) and (A16) 
require that 


g*|xi,V(xi),e( x ijj = g*|x 2 ,V(x 2 ),e(x 2 yj =0 (A22) 

Use of equations (A20), (A21), and (A2) shows that the relations in equation (A22) are 
identically satisfied. Equation (A14) now reduces to 


= 0 

90 dx^ayy 


(A23) 


The use of equations (A19) and (A20) in equation (A23) gives the differential 
equation 


v dV + eE = o ( A24 ) 

dx m 

which must be satisfied by V(x). Equation (A24) and the constraint equation 

N = 2A0^x 2 ) (A25) 

are sufficient to determine x^, x 2 , and V(x) where V^x-jJ = V^x 2 ) = 0. 

d f) 

By noting that 0 — 0 and V = — — 0' transform equation (A14) into the custo- 

dx 

mary Euler -Lagrange equation, the Legendre test (ref. 22) may now be applied to the 
second variation of the integral W*. The Legendre test states that in order for the 
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solution V(x) of equation (A24) to be such that the integral W* is maximized, the 
quantity 



2AmV 


must be negative definite over the interval from x^ to X 2 - This quantity is seen to be 
positive definite over the interval so that the extremum represented by V(x) cannot be 
a maximum, but must be either a minimum or a simple inflection. Other evidence pre- 
sented herein indicates that the extremum is a minimum. 
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Figure 3.- Variation of electron density for several values of a. 
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Figure 5.- Equilibrium waterbag contours for several values of a. 
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Figure 9.- Illustration of a multi p!e-contour-waterbag distribution. 
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(a) Time development from t = 0 to t = 1.25a^ 0 -1 . 

Figure 14.- Time development in phase space for a 1000-electron system near equilibrium. 
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(b) Time development from t = 2.50cc>p O " 1 to t = 5J5a)^ 0 ~ l . 
Figure 14.- Concluded. 
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